Intro to BioEnergeticFoodWebs

by Chris Griffiths, Eva Delmas and Andrew Beckerman, Dec. 2020.

1.8 ms

This document follows on from 'Julia in VS Code #1, #2 and #3' and assumes that you're still working in your active project.

This document introduces the BioEnergeticFoodWebs.jl and EcologicalNetworks.jl packages. It demonstrates how to run the BioEnergetic Food Web (BEFW) model, how to vary variables of interest (e.g., productivity) and construct experiments designed to investigate the effect of different variables on population and community dynamics.

For those that are unfamilar with the BEFW and it's application in Julia, we advise checking out the MEE paper before we start. Remember, the BEFW model is also based on a system of differential equations and is solved using the same engine as the DifferentialEquations.jl package.

10.5 μs

Packages

11.9 μs

First, import your package manager, activate and instantiate (only neccessary if you've closed VS Code between tutorials):

5.9 μs
5.3 s

Then install the BioEnergeticFoodWebs.jl, EcologicalNetworks.jl and JLD2.jl packages and let Julia know which packages you want to use:

6.8 μs
17.4 s

The JLD2.jl package will be useful later as it allows you to directly export and load a BEFW output object. Let's also set a random seed for reproducibility:

7.1 μs
MersenneTwisterseedstatevalsintsidxF
1002
idxI
0
19.2 μs

Preamble

One of main advantages of running food web models in Julia is that simulations are fast and can be readily stored in your active project. With this in mind, make a new folder in your project called out_objects (right click > New Folder). Alternatively, you can create an out_objects folder directly using mkdir().

8.0 μs

Running the BEFW

There are four major steps when running the BioEnergetic Food Web model in Julia:

  1. Generate an initial network

  2. Fix parameters

  3. Simulate

  4. Explore output and plot

Initial network

Before running the BEFW model, we have to construct an initial random network using the niche model. The network is characterised by the number of species in the network and its connectance value. Here, we generate a network of 20 species with a connectance value of 0.15:

11.9 μs
20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  0  0
 0  0  0  0  0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
 0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
32.8 ms

You can check the connectance of A using:

6.7 μs
co
0.145
2.9 μs

Parameters

Prior to running the BEFW model, you have to create a vector of model parameters using the model_parameters function. Numerous parameter values can be specified within the model_parameters function, however, most of them have default values that are built into the BioEnergeticFoodWebs.jl package. For simplicity, we use the default values here:

8.0 μs
p
Dict
1.0
:e_carnivore
0.85
:Γh
:m_producer
1.0
:c
0.0
:h
1.0
:vertebrates
:tmpA
:rewire_method
:none
:trophic_rank
:w
20×20 Array{Float64,2}:
 0.0  0.0  0.0   0.0   0.0   0.0   0.0       …  0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.0    0.0    0.0    0.0    0.0
 0.0  0.0  1.0   0.0   0.0   0.0   0.0          0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.0    0.0    0.0    0.0    0.0
 0.2  0.2  0.2   0.2   0.2   0.0   0.0       …  0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.25  0.25  0.25         0.0    0.0    0.0    0.0    0.0
 ⋮                           ⋮               ⋱  ⋮                           
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.125  0.125  0.125  0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.166667  …  0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.25         0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.0    0.0    0.0    0.0    0.0
 0.0  0.0  0.0   0.0   0.0   0.0   0.0          0.125  0.125  0.125  0.125  0.125
 0.0  0.0  0.25  0.25  0.25  0.25  0.0          0.0    0.0    0.0    0.0    0.0
:TSR_type
:no_response
:e_herbivore
0.45
:productivity
:species
:efficiency
20×20 Array{Float64,2}:
 0.0   0.0   0.0   0.0   0.0   0.0   0.0   …  0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.45  0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.45  0.45  0.45  0.85  0.45  0.0   0.0   …  0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.45  0.85  0.85     0.0   0.0   0.0   0.0   0.0   0.0
 ⋮                             ⋮           ⋱        ⋮                       
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.85  0.85  0.85  0.85  0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.85  …  0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.85     0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.0   0.0   0.0   0.0   0.0   0.0
 0.0   0.0   0.0   0.0   0.0   0.0   0.0      0.85  0.85  0.85  0.85  0.85  0.85
 0.0   0.0   0.45  0.85  0.45  0.85  0.0      0.0   0.0   0.0   0.0   0.0   0.0
:K
1.0
:S
20
:extinctionstime
:is_producer
:dp
#10
:x
:Z
1.0
:extinctions
:dry_mass_293
:np
5
:A
20×20 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 1  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0
 ⋮              ⋮              ⋮              ⋮           
 0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1  0  0
 0  0  0  0  0  0  1  1  1  1  1  1  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  1  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
 0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
:ar
:bodymass
:is_herbivore
more
942 μs

For more information and a full list of the parameters and their defaults values type ?model_parameters in the REPL.

6.8 μs

Simulate

To run the BEFW model, we first assign biomasses at random to each species and then simulate the biomass dynamics forward using the simulate function:

8.2 μs
Dict
:p
:B
8017×20 Adjoint{Float64,Array{Float64,2}}:
 0.524056  0.535102  0.977416     0.906369  …  0.459164  0.615476  0.245534  0.490372
 0.517543  0.527093  0.192875     0.923082     0.553908  0.718219  0.293288  0.625466
 0.50712   0.515335  0.0200148    0.710628     0.614887  0.701168  0.358189  0.746566
 0.501974  0.509129  0.00282734   0.491415     0.639262  0.621707  0.440936  0.839793
 0.505704  0.51205   0.000635531  0.328812     0.626891  0.530129  0.540591  0.893063
 0.518701  0.524421  0.000204033  0.217314  …  0.581946  0.441034  0.654733  0.901809
 0.539661  0.544868  8.4954e-5    0.144        0.51424   0.358944  0.77923   0.870207
 ⋮                                          ⋱                                
 0.559053  0.559053  1.60899e-5   0.0          0.0       0.139536  0.0       0.356004
 0.559054  0.559054  1.54769e-5   0.0          0.0       0.139536  0.0       0.356005
 0.559054  0.559054  1.45103e-5   0.0          0.0       0.139536  0.0       0.356006
 0.559055  0.559055  1.33312e-5   0.0          0.0       0.139536  0.0       0.356007
 0.559055  0.559055  1.21177e-5   0.0       …  0.0       0.139536  0.0       0.356008
 0.559055  0.559055  1.10852e-5   0.0          0.0       0.139536  0.0       0.356008
:t
375 ms

The simulate function requires the model parameters p and the species biomasses bm. In addition, you can specify the timespan of the simulation (using the start and stop arguments), fix a species extinction threshold (using extinction_threshold) and select a solver (using use). For more information type ?simulate in the REPL.

6.8 μs

Output and plot

Once the simulation finishes, the output is stored as a dictionary called out. Within out there are three entries:

  1. out[:p] - lists the parameters

  2. out[:B] - biomass of each species through time

  3. out[:t] - timesteps (these typically increase in 0.25 intervals)

The biomass dynamics of each species can then be plotted. Similar to the DifferentialEquations.jl package, the BioEnergeticFoodWebs.jl package also has it's own built in plotting recipe:

10.3 μs
22.7 ms

You'll notice that the biomass dynamics are noisey during the first few hundred time steps, these are the system's transient dynamics. The dynamics then settle into a steady state where the system can be assumed to be at equilbirum. You'll also notice that some species go extinct and some persist, the number of species in the food web can found using out[:p][:S] and the identity of those that went extinct using out[:p][:extinctions].

The BioEnergeticFoodWebs.jl package also has a range of built in functions that conveniently calculate some of the key metrics of the food web, these include the total biomass, the diversity, the species persistence and the temporal stability:

7.8 μs
biomass
2.0796559291208965
61.6 μs
diversity
0.7182050496471224
470 μs
persistence
0.55
100 μs
stability
-0.04410056440077254
93.8 μs

Each of these functions will output a single value. This value is the average over the last 1000 time steps. For more information, use ? to access the help files on each function in the REPL (e.g., ?species_persistence).

7.2 μs

Variables

Once you've got the BEFW model running, the next step is to vary a variable of interest and rerun. For example, we might be interested in what affect a small change in Z (consumer-resource body mass ratio) has on the estimated food web and its biomass dynamics. The default value for Z is 1.0, but what happens if we increase it to 10.0:

8.5 μs
866 ms

Similarly, what happens if we also increase the carrying capacity (K) of the resource from 1.0 (default) to 5.0:

6.4 μs
145 ms

As you've probably guessed, the main message here is that many variables can be changed in the BEFW model and it's super easy to do so. In the next step, we take this one step further.

6.6 μs

Experiments

The next step is to construct a computional experiment designed to investigate the effect of different variables on population and community dynamics. To do this we construct a gradient of variables as vectors and then simulate the BEFW model multiple times using a loop. To illustrate this, we're going to reproduce example 1 from Delmas et al. 2016. The aim of this example is to investigate the effect of increasing K on food web diversity. In addition, we're also going to allow α (interspecific competition relative to intraspecific competition) to vary and repeat the experiment 5 times with 5 different initial networks.

First, we define the experiment by creating vectors of our variables and fixing the number of repetitions:

11.3 μs
α
371 ns
k
18.1 μs
reps
5
46.0 ns

We then create a dataframe to store the outputs:

6.3 μs
df
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
32.0 ms

and construct a while loop to generate the 5 unique initial networks, each of which contains 20 species with a connectance value of 0.15:

6.7 μs
1.7 ms

We can then run the simulations by looping, using nested for loops, over the unique values of α and K, as well as the 5 unique initial networks. After each simulation we will save each output object to our active project as a JLD2 file and store any output metrics of interest in our dataframe:

7.3 μs
60.1 s

We can then explore the outputs and plot our results. Here, instead of using the built in plotting recipe, we will construct a plot that matches figure 1 in Delmas et al. 2016. Specifically, we will plot food web diversity (y-axis) as a function of K (x-axis) and α (colour):

8.5 μs
variablemeanminmedianmaxnuniquenmissingeltype
SymbolFloat64Float64Float64Int64Int64DataType
1
1.0
0.92
1.0
1.08
3
0
Any
2
:K
5.0
5.0
5.0
5.0
1
0
Any
3
:network
3.0
1.0
3.0
5.0
5
0
Any
4
:diversity
NaN
NaN
nothing
NaN
9
0
Any
5
:stability
-0.099386
-1.3342
0.0
-0.0
8
0
Any
6
:biomass
6.96612
2.99178
5.0
14.2182
9
0
Any
270 ms
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
1
0.92
5.0
1.0
NaN
-0.0
5.0
2
1.0
5.0
1.0
NaN
-0.0
5.0
3
1.08
5.0
1.0
NaN
-0.0
5.0
4
0.92
5.0
2.0
0.928022
-0.00466869
2.99184
5
1.0
5.0
2.0
0.948159
-0.00468039
2.99178
6
1.08
5.0
2.0
0.977442
-0.00468403
2.99182
16.8 μs
αKnetworkdiversitystabilitybiomass
AnyAnyAnyAnyAnyAny
1
0.92
5.0
4.0
NaN
-0.0
5.0
2
1.0
5.0
4.0
NaN
-0.0
5.0
3
1.08
5.0
4.0
0.413271
-1.3342
7.86188
4
0.92
5.0
5.0
NaN
-0.0
5.0
5
1.0
5.0
5.0
NaN
-0.0
5.0
6
1.08
5.0
5.0
1.0
-0.0
10.0
10.8 μs
pl
31.1 ms
shp
970 ns
ls
490 ns
clr
5.0 μs
lbl
852 ns
2.7 ms
367 μs

Finally, we can save our dataframe as a .csv file:

6.4 μs
"My_data.csv"
335 ms